function chkderiv(tchfe,tchfeMAP,tchfeERRCOV,sv)

[numvars,numcomp] = size(sv.mean);

opts = optimoptions("fminunc",FiniteDifferenceType="central");

%%% -- score function 
fprintf('Checking score output from emstep ... \n')
function [ll,gr] = chkscore(parms)
    vadist = vecparm(parms,numvars,numcomp);
    [~,ll,~,scr] = emstep(vadist,tchfe,tchfeMAP,tchfeERRCOV);
    ll = -sum(ll);
    gr = -sum(scr,1)';
end

sv1 = vecparm(sv,numvars,numcomp);
checkGradients(@chkscore,sv1,opts,Display="on");

%%% -- grad function 
fprintf('Checking gradient output from emstep ... \n')
function [ll,gr] = chkgrad(parms)
    vadist = vecparm_unc(parms,numvars,numcomp);
    [~,ll,gr] = emstep(vadist,tchfe,tchfeMAP,tchfeERRCOV);
    ll = -sum(ll);
    gr = -gr;
end

sv2 = vecparm_unc(sv,numvars,numcomp);
checkGradients(@chkgrad,sv2,opts,Display="on");


end